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Abstract 



> 

', Liischer's local bosonic algorithm for Monte Carlo simulations of quantum field 



theories with fermions is applied to the simulation of a possibly supersymmetric 
^ Yang-Mills theory with a Majorana fermion in the adjoint representation. Combined 

Q\ • with a correction step in a two-step polynomial approximation scheme, the obtained 



algorithm seems to be promising and could be competitive with more conventional 
algorithms based on discretized classical ("molecular dynamics") equations of mo- 
D ! tion. The application of the considered polynomial approximation scheme to opti- 



mized hopping parameter expansions is also discussed. 



1 Introduction 

Among quantum field theories the supersymmetric ones play a very distinguished role: they 
have much less free parameters than a general renormalizable quantum field theory with the 
same set of fields, and show remarkable non-renormalization and finiteness properties (for 
general references see, for instance, [^). Even more special are supersymmetric Yang Mills 
(SYM) theories. The Yang Mills theory with a Majorana fermion in the adjoint representation 
is automatically = 1 supersymmetric in the massless case, at least according to perturbation 
theory. The supersymmetry constraints in = 2 SYM theory can be exploited to determine, 
under some reasonable assumptions, its exact solution in the low energy limit 0. 

The non-perturbative properties of supesymmetric quantum gauge field theories enjoyed 
particular interest already in the 80's (see the review and references therein). After the 
recent beautiful exact results of Seiberg and Witten ^ there is a revived interest because 
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one can hope to understand better on these examples some general non-perturbative phenomana 
as confinement and chiral symmetry breaking. 

It is obvious that numerical simulations for studying non-perturbative properties of super- 
symmetric quantum field theories would be both interesting and desirable. In this way the 
assumptions beyond the exact solutions could, in principle, be checked and the results perhaps 
extended. Nevertheless, up to now it turned out to be impossible to find a lattice formulation 
with exact supersymmetry. The only possibility for reconciling lattice regularization with su- 
persymmetry seems to be to formulate a more general non-supersymmetric model and recover 
supersymmetry in the continuum limit as a result of fine-tuning the parameters. This approach 
was pioneered by Curci and Veneziano in the simple case of = 1 SYM, where there is 
only one parameter to tune (namely, the Majorana fermion mass). More general cases can be 
expected to be dealt with similarly, as shown, for instance, in the case of = 2 SYM in ref. 
0. This way of "embedding" supersymmetric theories in more general non-supersymmetric 
ones is natural since in Nature supersymmetry is also broken. Besides, as argued recently 
many of the "exotic" dynamical features of supersymmetric theories survive, if small symmetry 
breaking terms are added. 

Since supersymmetry connects bosons and fermions, a numerical simulation also involves 
the numerical approximation of fermionic Grassmann integrals. This is a notoriously diffi- 
cult problem, which makes numerical simulations of e. g. QCD with dynamical fermions quite 
hard. Nevertheless, there are several known fermion algorithms which do the job for QCD, 
for instance, the popular Hybrid Monte Carlo algorithm [||. This is not directly applicable in 
supersymmetric models where Majorana fermions are important. Nevertheless, some related 
algorithms based on discretized classical "molecular dynamics" equations of motions, as for 
instance the one in ref. 0, can work. However, recently Liischer proposed a conceptually 
rather different local bosonic algorithm |T^, which is an alternative for QCD [|r^, and can also 
be extended to supersymmetric cases. As recently pointed out by Borigi and de Forcrand ||12[| , 
this algorithm gains on attractivity, if it is combined with a "noisy correction" step as proposed 
in the early days of fermion algorithm developments |jl3|. This local bosonic formulation and 
its combination with a noisy correction step is the basis of the fermion algorithm investigated 
in the present paper. The formulation and tests were done in A^ = 1 SYM with SU(2) gauge 
group, but the methods used are obviously more general: they can be applied in many super- 
symmetric models with A^ = 1 and A^ > 1 supersymmetry and, of course, also in other models 
containing fermions and not related to supersymmetry. 

The plan of this paper is the following: In the next section the lattice formulation of 
models with Majorana fermions is discussed on the example of A^ = 1 SYM. In section |^ the 
optimized polynomial approximation scheme is introduced, which is then used in the local 
bosonic algorithms defined and tested in section Section || is devoted to the discussion of the 
applications of the optimized polynomial approximations to hopping parameter expansions. A 
possible fermion algorithm based on optimized numerical hopping parameter expansion is also 
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briefly discussed there. Finally, section |^ contains some concluding remarks. 



2 Majorana fermions on the lattice 

In order to define the path integral for a Yang Mills theory with Majorana fermions in the adjoint 
representation, let us first consider the same theory with Dirac fermions. The lattice action 
in case of Wilson's lattice fermion formulation has been considered in P]. If the Grassmanian 
fermion fields in the adjoint representation are denoted by ijj^ and ip^., with r being the adjoint 
representation index (r = 1, .., A*"^ — 1 for SU(A''c) ), then the fermionic part of the lattice action 

X fj.= l 

Here K is the hopping parameter, the irrelevant Wilson parameter removing the fermion dou- 
blers in the continuum limit is fixed to r = 1, and the matrix for the gauge- field link in the 
adjoint representation is defined as 

Vrs,., ^ Vrs,.,[U] = 2TTiUlTr.U.,Ts) = K:,.^ = • (2) 

The generators = |Ar satisfy the usual normalization Tr(Ar.As) = |. In case of SU(2) 
{Nc = 2) we have, of course, Tj. = with the isospin Pauli- matrices r^,. The normalization 
of the fermion fields in ([l|) is the usual one for numerical simulations. The full lattice action is 
the sum of the pure gauge part and fermionic part: 

S = S, + Sf. (3) 

The standard Wilson action for the SU(Ai'c) gauge field Sg is a sum over the plaquettes 

^. = /5E(l-l^ReTrf/,,) , (4) 

pi ^ / 

with the bare gauge coupling given by /? = 2Nc/g'^. 

In order to obtain the lattice formulation of a theory with Majorana fermions, let us intro- 
duce the Majorana field components 

v[/(i) = ^(^ + C^^) , ^^'^ ^ ^(-^ + C^^) (5) 
with the charge conjugation matrix C. These satisfy the Majorana condition 

^0) = ^0)T^ (j = l,2). (6) 

The inverse relation of (j^) is 

^ = ^(^(1) + , ip, = c^^ = ^(^(') - . (7) 

\/2 V 2 
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In terms of the two Majorana fields tlie fermion action Sj in eq. (|l]) can be written as 

^ X j=i ^,=l L J 

For later purposes it is convenient to introduce the fermion matrix 

4 

Qyv,xu — Qyv,xu[U] = Syx^vu ^ ^ ^ |^'^3/,2:+/i ( f ~l~ ^ fi)^vu,xfj, ~l~ ^y+fi,x(,^ 7/^)^ 



vu,yii 



(9) 



Here, as usual, fi denotes the unit vector in direction /i. In terms of Q we have 



^/ = E i^lQy^,xur. = ^ E E ^y^''Qyv,xu^^J^'' 
xu,yv ^ j=l xu,yv 

and the fermionic path integral can be written as 



(10) 



[d^d^]e-^^^ = detQ = n / W^^^]e 



(11) 



Note that for Majorana fields the path integral involves only [d^'^-'^] because of the Majorana 
condition in (|^). This implies for \[' = \['(^) or ^ = 



±VdetQ . 



(12) 



For a given gauge field the sign can be taken by definition to be positive, but for different gauge 
fields one has to care about continuity. The sign convention can be fixed, for instance, at the 
trivial gauge field Ux^j. = 1 to be positive. Then the positive sign stays by continuity until one 
reaches gauge fields with det(5[[/] = 0. 

The squre root of the determinant in eq. (|I2[) is a Pfaffian. This can be defined for a 
general complex antisymmetric matrix Mq,^ = —Mp^ with an even number of dimensions 
(1 < a,j3 < 2N) by a Grassmann integral as 

pf (M) ^ J [dc^]e-^-'''^^f^ = ^e,^^^...„,^,M„,^, . . . M^^p^ . (13) 

Here, of course, [dcf)] = d(j)2N ■ ■ ■ d(f)i, and e is the totally antisymmetric unit tensor. Using the 
above trick with doubling the number of Grassmann variables one can easily show that 

-I 2 

= det M . (14) 



pf(-M) 



m 



This establishes the connection between Pfaffians and determinants. According to eq, 
case of the above fermion action the antisymmetric matrix is Q = CQ. 

It is now clear that the fermion action for a Majorana fermion in the adjoint representation 

can be defined by 

4 

2 ^ 2 ^ 



X A'=l 



The path integral over ^ is defined by the Pfaffian pf(|CQ) = y'det(Cg) = ^det{Q). In 
this way the Majorana nature of \1/ imphes that one has to take the square root of the usual 
fermion determinant . Note that the Pfaffian defined by ([T3|) assigns a unique sign to the path 
integral for Majorana fermions, therefore there is no sign ambiguity. For positive determinant 
the Pfaffian is real, but if the determinant is negative then the Pfaffian has to be pure imaginary. 

Expectation values of Majorana fermion fields can be calculated as follows. For the Dirac 
fermion fields ip, ip we have, as is well known, 

i^v.^.^i^y.'^., ■ ■ ■ = I Mf/]e-^«[^] det Q[U] 

■ E 4\ll-ZQPV^Mu]^L ■ ■ ■Q[u]z\. > (16) 

with the antisymmetrizing unit tensor e and 

Z = J [dU]e-^^^^^ det Q[U] . (17) 

In order to express the expectation value of Majorana fermion fields by the matrix elements 
of the propagator Q[f/]~^, one can use again the doubling trick: one can identify, for instance, 
v|/ = v|/(i) and introduce another Majorana field in order to obtain a Dirac field. Then 
from eqs. (H) and (0) we obtain 

M/,,^,,vl/^^^,.^ . . . vl/^^^^,^) = Z,-/ J [dU]e-'^^^^^ det Q[U] 



■ (det Q[U])-' j [d^#]e-^«n^,, + CO(^,, + C^l^l) ■ ■ ■ (^,„ + C^Jj2-" , (18) 
where now 

Zm = l^f/le-^^l^VdetQif/] . (19) 
In particular, for n = 1 we have with Q = Q\U] 

(vp,^.) = Zll j [dU]e-'^^^^^ det Q[U] \ [Q'^I + C-^Ql^^C] . (20) 
The important case n = 2 can be expressed by six terms: 

■7 E {_^yiy2^ zi^l^ Z2X2 + ^xi ^x\y2Q ziyiQ Z2X2^yi + ^X2 ^y[x2Q z^x^Q Z2y2^y2 
Z\Z2 

^^xi Cx2 ^xiX2QziyiQz2y2^y\Cy2 '~ ^xi ^y\x^iQ Z1X2Q Z2y2^y2 ~ ^ X2 ^^2X2^ z^x^Q Z2yi^yi\ ' (^1) 

The indices on the charge conjugation matrix C show how the Dirac indices have to be con- 
tracted. 

Since the path integral over Majorana fermion fields is defined by a Pfaffian, the sign in eq. 



T2| ) is uniquely determined. It is easy to show that the fermion matrix in (|^) satisfies 

= , (22) 
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therefore det(Q) is always real. In principle, in numerical simulations one can take into account 
the phase of the Pfaffian by performing the Monte Carlo integration with the positive square 
root +^^1 det Ql, and correct for the phase by putting it into the expectation values. Although 
this solves the phase problem in principle, in practice there is still a problem because the 
numerical evaluation of the Pfaffian from the formula (|I3p is rather cumbersome. The only 
practical possibility in a Monte Carlo process is to monitor the lowest eigenvalues of Q and/or 
explicitly calculate det Q, in order to make it probable that det Q does not cross zero, where the 
question of the sign change in eq. (|T2p and the problem of the non-trivial phase of the Pfaffian 
emerge. In what follows it will be assumed that one can restrict the Monte Carlo simulation 
to gauge configurations where pf(|CQ) = +\J\ det{Q)\. 

3 Optimized polynomial approximations 

An important ingredient of Liischer's fermion algorithm [|T0| is the polynomial approximation 
of the function 1/x in some positive interval x G [e. A] (0 < e < A < oo). This can be 
achieved by Chebyshev polynomials, which minimize the maximum relative error ("infinity- 
norm"). Here we shall follow another approach by minimizing the deviation in an average 
sense. The problem can then be reduced to a simple quadratic minimization. The advantage 
of this approach is its fiexibility in choosing weight factors and regions of minimization. In 
addition, it is easily applicable to more general cases including broad classes of functions and 
also two-step approximations needed in the next section. In fact, the function 1/x refers to 
Nf = 2 degenerate quark fiavours. In the general case one has to approximate the function 
gince in the path integrals for one species of Majorana fermions the square root of 
the fermion determinant appears (see previous section), a Majorana fermion is obtained for 
Nf = 1/2, hence we have to approximate x~^^^. This corresponds to writing 

lVdet(g)| = {det(QtQ)}V4 ^ {det(Q2)}^/^ . (23) 
Here Q denotes the Hermitean fermion matrix defined by 

g = 75g = gt . (24) 

(See eq. (^.) 

In what follows, we shall always assume that in the numerical simulations the spectrum of 
g^g is bounded from below, and is contained in the interval [e. A]. In principle, this may induce 
some restrictions on the bare parameters (3, K, but in practical numerical simulations, e. g. in 
QCD, this condition seems always to be satisfied. Before discussing in section H the ways of 
implementing the fermion determinant (^) in the path integral, let us introduce some classes 
of the necessary optimized polynomial approximations. 
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3.1 Single step approximation 



In order to define an optimized polynomial approximation P{x) of the function = x^^-'/^ 
in the above interval, the first thing to do is to define a positive norm S characterizing the 
deviation of two functions from each other. The approximation is optimal if 6 is minimal. A 
simple possibility would be to minimize, for instance. 



dxx^ 



X-" -P(x] 



(25) 



Here x'^ is some arbitrary power giving different weights to different parts of the interval. 

In what follows, we shall only consider the special case w = 2a, which corresponds to 
considering the relative deviation. This gives the deviation norm 



(26) 



The square root in the definition of the norm is a matter of convention. One can, as well, 
minimize 



52 = £ dx [1 - x"P(x)]^ = ^ 



.l+a+n—u A+ct+n—v n n 



X 



dx 



u=0 



1 + a + n — u 



ui=0 U2=0 



1 + 2a + 2n — Vi — V2 



■ (27) 



The coefficients of the polynomial Cy are assumed to be real. 

The second form in ( pTj) shows the advantage of choosing norms of the type (p5|)-p^): since 
52 is a quadratic form of the coefficients, the minimum can be easily determined to be at 



1/1=0 



(28) 



where from 



we have 



X 



1+a+n—v 1+a+n—v 



{n)u 



M, 



{n)v2Ui 



M, 



1 + a + n — u 

y^l+2a+2n—u\—U2 ^l+2a+2n—v\—y2 



(n)u\U2 



1 + 2a + 2n — ui — U2 
The optimized polynomial approximation is then 

n 

Pnia; e, X;x) = J2 ^nui^; e, A)x""'' . 

i/=0 

The roots of this polynomial will be needed in the form 

r„j(a; e. A) = rj = {fXj + ivjf = ^] - v] + 2%[ijVj (j = 1, . . . , n) 



(29) 



(30) 



(31) 



where Vj > will be assumed. Since the coefficients Cy are real, the roots come in pairs with 
equal and opposite fij. One can then write 



P„(a;e, A;x) = c„o(«; e. A) J| [x - r„j(a; e. A)] = CnoY[[iVx + Hjf + u"^] . 



(32) 



For a good approximation usually quite high orders n ^ 1 are necessary. In this case 
the numerical determination of the coefficients Cnuicn; e, X) and roots rnj{(y;e,\) can become 
cumbersome. In particular, very high numerical precision is necessary, which is not available 
e. g. in Fortran. However, algebraic programs usually provide the possibility to calculate with 
arbitrarily high precision. For instance, using this feature of Maple, one can obtain the required 
precision, and use the results as input data for Fortran programs. For the calculation of the 
inverse matrix in (|28|) one can use LU-decomposition with pivoting |jl4|. The roots of the 
polynomial can be numerically determined by Laguerre iteration . 

For gluinos on the lattice we need the power a = 1/4. In this case good approximations can 
be achieved by polynomials of order n = 32 or n = 40. An example for the spectral interval 
[e = 0.03, A = 4.0], which is typical in the test runs discussed in section ^^3] , is illustrated in 
figure m. 



relative deviation from x'^ (-1/4) : n=40 



. 0001 - - 




^2 , 5 




-0 .0001 



-0 . 0002 



-0 .0003 



Figure 1: The relative deviation of P4o(^; 0.03, 4.0; x) from x The value of the 
deviation norm is here: 6 = 2.51.. ■ 10"^. 
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3.2 Two-step approximations 

In the next section two-step approximations will be used for the local bosonic description of 
the fermion determinant. Let us assume that a first polynomial approximation P{x) of the 
function is already known. This can be obtained, for instance, by the method described 
in the previous subsection for some relatively low order n. The task is now to find the best 
approximation in the product form P{x)P{x), with a polynomial P{x) of order n > n. 

This can be done, in fact, quite similarly as in subsection lO. The equation replacing eq. 



(|27|) is now 

6^ = f dx 



1 - X! H CvCuX^ 



a+n+n—v—y 
j/=0 p=0 J 

Here Cp denotes the coefficients in P. Everything remains the same as in subsection |3.1|, only 



(33) 



(|29|) has to be replaced by 

1+a+n+n—u—u 

n n y^+2a+2n+2n—v\—V2—v\~V2 ^l+2a+2n+2n—u\—V2—v\—U2 

The coefficients of the optimized polynomial P{x) = Pn,n{ci; e. A; x) are given by 

n 

Cn,nu{0i; e, \) = ^{n,n)uui^{n,n)u, ■ (35) 

A somewhat different type of approximation will also be needed when, for the given poly- 
nomial P{x), the optimized P{x) is searched which is minimizing 



6 = dxx 
In this case, instead of (|3^), we have 



,A _ 2 



(36) 



n ^l+w+a+n+n—u—u ^l+w+a+n+n—i'—i' 

V(n,n)u y, Cp - ■ ■ ■ ■ - - , 

1 + w + a + n + n — u — u 

^l+w+2n—ui—U2 gl+«i+2n— i^i — 

M(^n,n)u2Ui = 7—. "7^ ■ (37) 

1 + w + 2n — ui — 1^2 
3.3 Approximations in the complex plane 

The above polynomial approximation scheme can also be extended to the complex plane. De- 
noting the complex variable by {x + iy), it is possible to approximate the function {x + iy)~°' 
in some region of the complex plane. This problem occurs, for instance, in the non-hermitian 
variants of local bosonic algorithms |jl2[. The approximation regions have to cover the complex 



spectrum of the fermion matrix Q. Since, according to the relation (p2|), the eigenvalues of Q 



come in complex conjugate pairs, the region of approximation is symmetric with respect to a 
reflection on the real axis. Natural shapes of the complex region are circle or ellipse, but a 
simple symmetric rectangle is also appropriate. 

In this latter case, for instance, one can minimize 



dx / dy\x + iy[" {x + iy) ^ — P{x + iy) 



(38) 



Therefore, for w = 2a, eq. ( p7D is replaced by 



dx I dy 



^c^{x + iy) 



a+n—u 



v=0 



(39) 



In the interesting cases, when a = Nf = |, 1, 2, . . ., the integrals can be performed analytically 
and give rise to relatively simple expressions in A,e,7. The coefficients c„y(a;e, A,7) of the 
optimal polynomial approximation P„(q;; e. A, 7; x + iy) are given again by a similar expression 
as eq. (pSf) . 

An illustration for the positions of the roots of the polynomial P„ for a = 1, (e = 0.1, A = 
2.0, 7 = 1.0) is figure |^. The spectral region is typical for the situation in the tests discussed 
in section [4.3| , similarly to fig. |l|. As one can see, the roots wind around the rectangular region 
in the complex plane. The value of the deviation norm is 5 = 0.0132... This is relatively 
large for this high order, showing that the rectangular region is not optimal for this kind of 
approximations. Elliptical shapes are more natural from the mathematical point of view. A 
rectangle is, however, more convenient from the point of view of monitoring the eigenvalue 
ranges. For instance, one can easily determine the extremal eigenvalues of the Hermitean 
matrices {Q + Q^) and i{Q^ — Q), and infer from them the necessary rectangle. (The extremal 
eigenvalues can be determined, for instance, by gradient methods ||l6l.) In addition, as a closer 
inspection shows, most of the contribution to 6 comes from the areas of the rectangle near 
the corners, where there are typically no eigenvalues anyway. In fact, in the important central 
regions the approximation is several orders of magnitude better. 



4 Local bosonic algorithms 

Let us now turn to the question, how to implement the fermion determinant factor in eq. (|23| ) 
in the path integrals appearing, for instance, in the expectation values (p!8|)-(|21|). Following 
ref. [0 and taking some polynomial approximations, as the Chebyshev polynomials or the 



polynomials introduced in section |3.1| , we have 

{detg2}i/4 L_ = L_ . (40) 

det[P„(g2)] det(c„o) n"=idet[(g + /i,)2 + i/2] 

Here — ^ means a polynomial approximation with deviation norm equal to 6. The determinant 
factors can be written with the help of complex pseudofermion fields (pjx, {j = l,..,n) as 
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polynomial roots for x'^(-l) : n=40 



1-- 



Y 

0.5^ 



0.5 



1.5 



Figure 2: The positions of the roots of P4o(l; 0.1, 2.0, 1.0; x + iy) in the complex plane. 



Gaussian integrals: 



det[(g + /ij)2 + 



oc / exp - 5] (l)*y[iy]6y, + iQ + 



(41) 



xy 



This is how the fermion determinant is represented in path integrals by a local bosonic action. 



4.1 Two-step approximation with noisy correction 

In order to obtain a very high precision approximation one has to take a large number of 
pseudofermion fields. This increases the required storage space and might also cause problems 
with long autocorrelations. The solution is to introduce a two-step approximation scheme, as 



discussed in section 3.2 



In the first step we make an approximation by the polynomial P{x) with a deviation norm 



6: 



{det Q'} 



211/4 _J_ 



det [P(g2)] 



(42) 
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This is improved to a smaller deviation norm 5 in the second step by multiplying with the 
polynomial P{x): 

{detQ2}i/4 ^ _ . (43) 

^ det[P(Q2)]det[P(Q2)] ^ > 

The first approximation is realized by pseudofermion fields as in eq. (0), whereas the correction 
factor det[P((5^)] is taken into account by a Metropolis correction step, similarly to ref. |jl3 



Since this latter is realized by Gaussian i. e. "noisy" unbiased estimators, we shall call it here 
briefiy noisy correction. 

By the Monte Carlo updating process one has to reproduce the canonical gauge field distri- 
bution 

w[U] = ^ = , (44) 

^ ^ det[P(Q[f/]2)]det[P(Q[?7]2)] ^ ^ 

where Sg[U] is the pure gauge field part of the action in Using ( ^T]) to represent det P, the 
canonical distribution in terms of the gauge field U and pseudofermion field can be written 
in the form 

w[U, 0] = = . (45) 

^ det[P(Q[f/]2)] ^ ^ 

The updating of the field can be straightforwardly done by heatbath and/or overrelaxation 
algorithms. 

In order to prove detailed balance for the updating of the gauge field, we have to show that 
the transition probability P{[U'] [U]) of the Markov process satisfies 

Pm - [U']) ^ e-^.^[^]det[P(g[[/f)] 

P{[U']^[U]) e-V[f^']det[P(Q[[/]2)] ■ ^ ^ 

The transition probability is given by 

P{[U'] - [U]) = Pm{[U'] ^ [U])PA[U'] ^ m , (47) 

where Pm{[U'] ^ [U]) is the transition probability of the updating of the gauge field with 
the action Sgp[U], and Pa{[U'] ^ [U]) is the acceptance probability of the correction step. 
P]\i{[U'] <— [U]) is realized by standard Metropolis sweeps. It is assumed to satisfy detailed 
balance for the action Sgp[U], that is 

PmHU'] ^ [U]) e-^.^[^'l ■ ^ ^ 

This can be achieved in different ways, for instance, by a randomly chosen order of links or by 
updating always one half of the links in a checkboard decomposition or by doing the sweeps 
pairwise with opposite orders of link updates. Comparing eqs. (^61) -(^HD one can see that the 
condition for the acceptance probability is 

PA[U'] ^ [U]) det[P{Q[Un 



PaHU] - m) det[P(g[t/f )] 
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(49) 



This could be satisfied by the the standard Metropohs acceptance probabihty 

P.(if/']^[f/])^nnn|l,i»^|. (50) 
J ^ \ 'det[P(Q[f/']2)] J ^ ^ 

However, this would be too expensive to implement numerically. 

The idea of the noisy correction is to generate a random vector t] according to the 
normalized Gaussian distribution 

g-7?tp(Q[t/]2)r, 

J[dr]]e-^'P(Ql^^')v ' ^^^^ 
and to accept the change [U'] ^ [U] with probability 

min{l,A(r/;[f/']^ [[/])} , (52) 

where 

A{r^; [W] ^ [U]) = exp [-r^^ {P(Q[Uf) - P{Q[Ur)}v) • (53) 
This means that the acceptance probability is 

Due to the relation 

the detailed balance condition in ( ^9|) is satisfied. 

In this form the algorithm would be still too expensive, because of the necessity to generate 
the distribution in eq. (|5TD . One can, however, easily generate the simple Gaussian distribution 



(56) 



and then obtain the above r] from 

= P{Q{Uf)-^r,' . (57) 

The exact calculation is still too difficult, but one can obtain the result to a required approx- 
imation by the two-step polynomial approximation technique introduced at the end of section 
^ An approximation with the same deviation norm 5 as in (|43|) is sufficient. 

According to eq. (|i3| ) the polynomial approximation R{x) for P{x)~^^'^ has to obey 

R{x) ^ P(x)-5 x5P(x)^ ^ x^^S{P{x)) . (58) 



Here we defined the approximation 



S{P) ^ P^ , (59) 
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which can be obained by a single-step polynomial approximation defined in section |3.1| for 
a = — I in the interval [A"^/^, e"^/^]. After this, R{x) in can be determined by minimizing 
the deviation norm 6 defined in eq. (|36D for a = 1/8 and P{x) replaced by S{P{x)). The weight 
factor is somewhat arbitrary, but a good choice is, for instance, w = 0. This means that in the 
interval [e. A] the absolute deviation is minimized. 

This algorithm based on two-step approximations is not "exact" because for any finite S 
there is still some deviation from the functions required in eqs. ( ^31) and (^8|-^). This is 
in contrast to the exact algorithms with noisy correction proposed for QCD in refs. [O, ITB 



It is, however, not clear how to extend the exact algorithms to the case with square roots 
considered in this paper. Nevertheless the two-step approximations can easily be made very 
precise. Namely, provided that the second approximation is already good, the acceptance rate 
in the noisy correction step mainly depends on the first polynomial in (|4^), and only very 
little on the precision of the approximation achieved by the second one. The application of the 
polynomials in eq. (^) and (^) is not very time consuming, therefore the decrease of 6 does 
not lead to a dangerous slowing down of the simulation. This applies both to the gluino model 
considered in this paper and to QCD with quarks in the fundamental representation. Therefore 
the two-step algorithm proposed here is also a potential candidate for simulating QCD. 



4.2 Noisy estimators 

The pseudofermion fields can also be used to obtain the matrix elements of the fermion prop- 
agators. The required expectation values are easily obtained during the updating process, and 
in this way the fermion matrix elements can be determined. The expressions in terms of the 
pseudofermion fields are usually called noisy estimators. 
The starting relation is, for a general matrix M, 

det M = M^-J det M . (60) 

In our case we have from this 



-[detQ 



211/4 



[detg^]V4 • ^^^^ 



Combining this formula with eqs. (|40|)-(^TD we obtain 



i=i 

Here the notation 

= Y.(Q + ^^j)^y(f^jy (63) 
y 

is used. For instance, for the gluino condensate eqs. (pO|) and ( |62D give 

n 

- (i^.^.)) = -2 ^ ((0*.750, J + (0;.750,.))^^ . (64) 
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(Note that this is the unrenormalized condensate, which has to be subtracted and renormahzed, 
in order to obtain a physical quantity.) 

The expectation values of several pairs of fermion variables can also be obtained in a similar 
way. For instance, for two pairs one can use 

dM^pdM^s 

and in our case 

4/92 &' 



det M = (m^-JM,-i - M-^^M^^) det M , (65) 



) [det Q'Y" = {Q-^'Q:^ - Q;J,Q-^) [det QT^' • (66) 



3 \dC^yxdC^ xvz ^Qwx^Qyz / 

As in (^), it is straightforward to transform this in an expectation value in terms of the 
pseudofermion field. 



Table 1: The values of some simple characteristic quantities on 4^ ■ 8 lattice 
at (3 = 2.0, K = 0.150. The quantities are defined in the text. In the sixth 
column the number of combined sweeps is given in thousands. For the two-step 
algorithms the seventh column shows the updating probability of gauge links, 
and the last column contains the acceptance rate of the noisy correction step. 



{n, )n 


□ 


\Pi\ 


-(^.^.) 


X2 


sweep 


Pu 


acc. 


2 


0.50068(19) 


0.05088(23) 


9.4887(24) 


0.54831(93) 


28 






4 


0.51117(16) 


0.05095(23) 


11.1486(45) 


0.5668(10) 


40 






8 


0.50857(38) 


0.05191(51) 


11.818(12) 


0.5698(24) 


20 






16 


0.50675(8) 


0.0488(8) 


11.908(24) 


0.5657(39) 


20 






24 


0.5075(14) 


0.0524(16) 


11.897(45) 


0.5619(72) 


15 






4,24 


0.50611(29) 


0.0506(3) 


11.160(4) 


0.5561(13) 


60 


0.1 


0.59 


6,24 


0.50578(32) 


0.0514(3) 


11.639(8) 


0.5563(14) 


25 


1.0 


0.54 


8,24 


0.50600(40) 


0.0508(6) 


11.833(12) 


0.5556(21) 


20 


1.0 


0.74 



4.3 Numerical simulation tests 

The first tests of the algorithm described in the previous subsections were performed in the 
SU(2) Yang Mills model with gluinos. The bare parameters were j3 = 2.0 and 0.125 < K < 
0.175. The lattice sizes were in the range 4^ ■ 8 to 8'^ ■ 16. 

In order to observe the n-dependence of some simple global quantities and the corresponding 
autocorrelations, a series of tests were done on 4'^ ■ 8 lattice at {f3 = 2.0, K = 0.150). This is 
in a reasonable range of parameters for numerical simulations, as the obtained estimates of 
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the value of the string tension in lattice units show (see below). For the spectral interval 
[e = 0.03, A = 4.0] was taken, which contains the eigenvalues of confortably: it turned out 
that in equilibrium configurations the eigenvalues were always larger than 0.05 and smaller than 
3.7 (they were monitored at every sweep by a gradient algorithm |16|). The chosen quantities 
were: the plaquette □ = W\;y = ^Trf/pj, the absolute value of the Polyakov line in the time 
{Lt = 8) direction |P/|, the unrenormalized gluino condensate from the noisy estimator in eq. 
(|6^) , and the k®k Creutz ratio with k = 2 calculated from the Wilson loops as 

yyk,k-iyvk-i,k 
Here a^cx is the string tension in lattice units. 

The results with 2 < n < 24 for the single-step algorithm, and with n = 4, 6, 8; n = 24 
for the two-step algorithm discussed in subsection [4.1| are shown in table |^. The corresponding 
integrated autocorrelations are given in table |^. 



Table 2: The integrated autocorrelations of different quantities expressed in 
numbers of combined sweeps on 4^ ■ 8 lattice at/3 = 2.0, K = 0.150. The 
quantities are defined in the text. 



{n, )n 


□ 


\Pi\ 


-(^.^.) 


X2 


2 


6.4(3) 


1.0(1) 


1.0(1) 


2.8(1) 


4 


13.9(3) 


1.4(1) 


1.5(1) 


5.2(1) 


8 


36(2) 


3.3(3) 


1.6(1) 


15(1) 


16 


117(14) 


8.3(3) 


1.9(1) 


35(4) 


24 


165(35) 


19(2) 


2.3(2) 


43(5) 


4,24 


29(3) 


4.0(3) 


1.6(3) 


12.8(6) 


6,24 


18(2) 


1.7(4) 


1.8(3) 


7.8(7) 


8,24 


23(4) 


2.7(5) 


1.7(4) 


9.1(6) 



No particular effort was made to optimize the mixture of the different updating steps. This 
should certainly be an important part of the final optimization [|1^, which should be done 
in larger scale applications. Here only the n-dependence was investigated for a fixed reasonably 
looking mixture: 1 heatbath and 6 overrelaxation sweeps for the pseudofermions followed by 2 
Metropolis sweeps with 8 hits per link for the gauge field. In case of the two-step algorithm 
before the accept-reject correction step 2 Metropolis sweeps were done with opposite orders of 
links. This was repeated 5-times. In order to reach a good acceptance, at every link it was 
first decided whether to update it or not. By choosing the updating probability < 1, the 
acceptance rate could be tuned appropriately. In table the numbers of "combined sweeps" 
always mean the numbers of these sequences of sweeps. 
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Tables |l| and ^ show that for the given statistical errors n = 24 is enough, and from the point 
of view of autocorrelations the algorithm with two-step approximation and noisy correction is 
better. In addition, compared to a single-step algorithm with the same approximation order n, 
the combined sweeps of the two-step algorithm also need less computer time. In table |^ it is 
amazing to observe that the physically more interesting Creutz ratio X2 has definitely smaller 
integrated autocorrelations than the plaquette □. The noisy estimators for — always 
have very short autocorrelations. 

In the two-step algorithm it is an interesting question, how the acceptance probability of 
the noisy correction is behaving for increasing lattice sizes. In order to see this, the simulation 
with n = 4, n = 24 was scaled up to 6^ ■ 12 and 8^ ■ 16, always with n = 24 in the correction 
step. It turned out that acceptances similar to n = 4 on 4^ ■ 8 could be reached with n = 6 and 
n = 8 on 6^ ■ 12 and 8^ ■ 16, respectively. With the same mixture of sweeps and link updating 
probability pu = 0.1, as in table |l|, the acceptance rates were 55% and 65%, respectively. This 
is good news, which tells that the number of pseudofermion fields for the first approximation 
has to grow at most linearly with the lattice extension. 

Another important question is the behaviour of these local bosonic algorithms for decreasing 
lattice spacing and/or fermion mass. A detailed investigation of this could not be carried out 
here. From the data in table [l| follows that the square root of the string tension in lattice units 
is at {(3 = 2.0, K = 0.150) roughly ay/a ^ 0.75. Defining the physical scale by y/a ~ 0.45 GeV, 
we obtain a ~ 1.7 GeV~^. In order to see the behaviour of the algorithm and of a^/a as a 
function of the fermion mass (K), numerical simulations on 6^ ■ 12 lattices were performed at 
(3 = 2.0 and K = 0.125, 0.150, 0.175 with the single-step algorithm. The approximations 
were done by the polynomials Ps{l/4:;e = 0.07, A = 3.5; x), Ps{l/4:;e = 0.03, A = 4.0; x) 
and P8(l/4; e = 0.005, A = 5.0; x), respectively. In all three cases the chosen spectral intervals 
comfortably covered the observed extrema of the spectrum of Q^. About 1000 combined sweeps 
were done after equilibration. The plaquette autocorrelations were about 10, 30 and 50 complete 
sweeps, respectively. At K = 0.125 and K = 0.150 the square root of the string tension in 
lattice units were roughly the same as the above value. At K = 0.175 a small decrease to a 
value ay/a ~ 0.7 was observed. 

Another question is the size of the effects of dynamical gluinos on the expectation values. 
In the considered range of bare parameters this is still small. A short simulation in pure gauge 
theory gives, for instance, at /3 = 2.0 on 8^ • 16 lattice □ = 0.50113(6) and X2 = 0.5490(2). 
These are only slightly below the values in table |l|. Larger effects are expected at larger j3 in 
the critical region of the hopping parameter K. 

The experience with the local bosonic algoritms on these relatively small lattices at rela- 
tively low lattice cut-offs was quite positive. As remarked above, from the point of view of 
autocorrelations and computational speed the two-step algorithm with noisy correction was 
clearly better. Of course, detailed tests on larger lattices and at smaller lattice spacings are 
necessary. A direct comparison with the conventional algorithms based on discretized classical 
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equations of motion would also be interesting. 



5 Other applications of the optimized polynomials 
5.1 Optimized solvers for matrix inversion 

The optimized polynomial approximations defined in section R| can also be used for matrix in- 



version instead of the Chebyshev polynomials in Chebyshev iteration schemes IjT^ • In numerical 
simulations of fermionic field theories an important task is to calculate the fermion propagator 
matrix elements in a given gauge field background (^[f/]^^ = Q~^. An approximate solution of 
the equation 

Qp = v (68) 

is the vector 

Pn = Pn{l;e,X,-f;Q)v = Cno(l;e,A,7) - r„j(l; e, A, 7)] j v . (69) 

Here we used the first form of the optimized polynomial P„ in (|32|), but other forms are also 
possible. Since we are actually considering the non-Hermitean fermion matrix Q, the approx- 
imation has to be optimized in the region of the complex plane x E [e, A], ?/ G [—7, 7], which 
covers the spectrum of Q. The deviation norm 6, which gives the precision of the approximation 
in (|69D, can be characterized by an expression similar to ( ^91) weighted by the spectral density 
of the matrix Q. 

For illustrative purposes figure ^ shows the squared length of the residue vector 

rn = V - QPn (70) 

obtained from (p9D on a typical 4^-8 configuration at {/3 = 2.0, K = 0.150). For comparison, the 
same quantity is also shown as a function of the number of iterations in the popular conjugate 
gradient (CG) algorithm. As the figure shows, an application of the optimized solver P32 reduces 
the residuum roughly as much as 32 CG iterations. 

An iterative scheme is obtained if the optimized polynomials are successively applied to the 
residuum vectors. The resulting cyclic iteration process can be represented by the following 
simple iteration equations: start by defining p^^ = p„ and r!^^ = rn and then set for /c = 1, 2, . . . 

r(^) = V - Qpi'^ , 
p^^ tends to the solution p of ( |68D in the limit k 00 because 

pi^') = PnV + (1 - PnQ)pL'-'^ = Pn EV " PnQY V • (72) 

i=o 
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Residuum in CG iterations 




10 20 30 40 

iteration 

Figure 3: The square of the residuum vector in a CG iteration on a 4^ ■ 8 configuration at 
{(3 = 2.0, K = 0.150). Twelve initial vectors are taken, starting from a randomly chosen 
point. The two upper bunches of curves belong to the normal CG iteration, whereas in 
the lower bunch the CG iteration is started after the multiplication of the twelve initial 
vectors by the optimized solver ^32(1; 0.1, 2.0, 1.0; x + iy) having 6 = 0.024... 

Therefore the iteration is realizing the simple identity 

= PniQPn)-' = Pn[l - (l " ^.Q)]"' • (73) 

For an illustration of this cyclic optimized solver (COS) iteration scheme see figure ^ where 
instead of Q the positive matrix Q^Q is inverted. In the figure the square of the residuum vector 
is shown as a function of the number of matrix multiplications (in this case a CG iteration step 
requires only one matrix multiplication). As the figure shows, the COS iteration with a high 
enough order polynomial has a similar performance as CG iteration. 

The application of the optimized solver P„ requires n matrix multiplications and no vector 
norm calculations. Therefore the amount of arithmetics in COS iterations compares favourably 
with the CG iteration. Nevertheless there are also known methods which are better than CG 
iteration [|I8] , therefore the best choice of the inversion method at given bare parameters and 
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Residuum in COS iterations 
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Figure 4: The square of the residuum vector in the cychc optimized solver iteration 
scheme for different polynomial orders n = 40, 60, 100. The vector v is always the same 
and the configuration is at /3 = 2.0, K = 0.175 on 8'^- 16 lattice. The interval of polynomial 
approximations is [0.01,4.4]. The dotted curve shows the result of conjugate gradient 
iteration. 

lattice sizes has to be decided on a case by case basis. 

The optimized solver can also be considered as an optimized hopping parameter expansion. 
Let us write the fermion matrix in (P) as 



Q=l- KM 
with the hopping matrix M. Then we have 



(74) 



P„(l; e. A, 7; Q) = c„o(l; e. A, i)Y[[l- r„, (1; e. A, 7) -KM] = Y. ^-^{1; e. A, ^)K''M^ . (75) 



i/=0 



This can be compared to the n'th order hopping parameter expansion of the inverse matrix 

n n 

J2 K^M" = (1 - K"+iM"+i) /(I - KM) = (-1)" n [e2^"/("+^) - KM] . (76) 

u=0 j=l 



20 



The only change in (|75D with respect to this is the different (optimized) choice of the roots of 
the n'th order polynomial. From the point of view of the numerical iterative hopping parameter 
expansion |]19|, this is practically no change in the arithmetics, but a considerable gain in 
precision. In addition, there is no convergence problem because for n — oo we always have 
5^0. 



5.2 An algorithm with optimized hopping parameter expansion 

The optimized hopping parameter expansion of the previous subsection can also be used in the 
old fermion algorithm based on the direct evaluation of the change of the fermion determinant 



2T| , pOf . Let us remember that according to the change of the fermion matrix is 

^Qyv,xu ^ {^y,z+p^zx(^'^ ~l~ 'yp)^^vu,zp ~l~ ^yz^z+p,x(^^ Tp)^^^f,zp] 5 ('''''') 

if the link zp is changed according to 

K'p = + A\4, . (78) 

The ratio of the two fermion determinants is 

det(Q + Ag) 



det(g) 



det(l + AQ ■ Q-') . (79) 



The 24 ® 24 matrix needed for the calculation of the determinant on the right hand side can be 
approximately calculated by the iterative numerical hopping parameter expansion, or by the 
improved variant discussed above using optimized solvers for the inverse. 



6 Discussion 

The local bosonic algorithms for gluinos on the lattice investigated in section ^ seem suitable 
for numerical experiments aiming at an understanding of Yang Mills theories with massive 
gluinos. The study of the limit when the gluino mass goes to zero should be possible and could 
reveal the character of the = 1 supersymmetric limit by comparing to the expectations based 
on previous analytical work In particular, the combination of Liischer's local bosonic 

approach with the noisy correction method |13] turned out to be rather effective in the 



performed tests. Since this combined algorithm has several algorithmic parameters, the final 
optimization for given lattice sizes and bare parameter values is a non-trivial task, which has 
to be carried out in the particular applications. 

The required polynomial approximations are defined in a scheme which is different from 
the usual approach leading to Chebyshev polynomials. The deviation norm is defined as an 
integral which is quadratic in the polynomial coefficients (see section |]). This allows a great 
flexibility in choosing weight factors and approximation regions. The optimal approximation 
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in this average sense is well suited for the present purposes, in particular, in the two-step 
approximation scheme applied in section ^. 

Besides the local bosonic algorithms, other applications of the optimized polynomial ap- 
proximation technique of section ^ have also been briefly discussed. In particular, the use of 
optimized solvers for matrix inversion and the application of the optimized hopping parameter 
expansion for a fermion algorithm without pseudofermion fields have been considered in section 
^ The test of this latter algorithm goes beyond the scope of the present paper. 

An important question in these algorithms is the role of the deviation norm 6, which is 
characterizing the quality of the polynomial approximations. Since exact results are obtained 
only in the limit 6^0, this could necessitate an extrapolation of the expectation values to S = 
0. Nevertheless, in the two-step approximation scheme a relative precision with 6 ~ 10~^ — 10~^ 
seems possible with reasonable effort, therefore the practical need for an extrapolation is not 
clear. For instance, such a precision comes already close to the machine precision on 32-bit 
computers, which is generally considered to be enough in present numerical simulations. 

Although the methods used here have been formulated and tested in case of the SU(2) 
gauge theory with a Majorana fermion in the adjoint representation, they are also applicable 
in other quantum field theories with fermion fields, as QCD or Higgs- Yukawa models including 
also scalar fields. 
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